home *** CD-ROM | disk | FTP | other *** search
/ NetNews Offline 2 / NetNews Offline Volume 2.iso / news / comp / lang / c-part1 / 6762 < prev    next >
Encoding:
Internet Message Format  |  1996-08-05  |  4.6 KB
  1. Path: keats.ugrad.cs.ubc.ca!not-for-mail
  2. From: c2a192@ugrad.cs.ubc.ca (Kazimir Kylheku)
  3. Newsgroups: comp.lang.c
  4. Subject: Re: Best way to make a circle
  5. Date: 14 Feb 1996 14:23:09 -0800
  6. Organization: Computer Science, University of B.C., Vancouver, B.C., Canada
  7. Message-ID: <4ftncdINNr8p@keats.ugrad.cs.ubc.ca>
  8. References: <1996Feb8.043040.19301@lafn.org> <4fo6pt$39o@worf.qntm.com>
  9. NNTP-Posting-Host: keats.ugrad.cs.ubc.ca
  10.  
  11. In article <4fo6pt$39o@worf.qntm.com>,  <pb@netcom.com> wrote:
  12. >In <1996Feb8.043040.19301@lafn.org>, an234@lafn.org (Andres Lessing) writes:
  13. >>
  14. >>I am working on 3d Graphics quite succesfully as of now and am trying to 
  15. >>put together a good Graphics library that is much faster than BGI drivers.
  16. >>I know that Sin and Cosine take to much time to calculate... so... which 
  17. >>way should I do it?
  18. >>
  19. >First a little math...
  20. >Circle: x = sin(t)
  21. >        y = cos(t)  step t from 0..2pi
  22. >
  23. >difer:  dx = +cos(t)/dt
  24. >        dy = -sin(t)/dt
  25. >
  26. >subst:  dx = +y/dt
  27. >        dy = -x/dt
  28. >
  29. >let dt=1/64 and we have in c
  30. >
  31. >        x += y/64;
  32. >        y -= x/64; /* rad/64 or about 400 steps around */
  33. >
  34. >if we let m = 2pi/steps we get...
  35. >
  36. >        x += m*y;
  37. >        y -= m*x;  /* 'tiz better to use the new x in this step */
  38. >
  39. >Using the above you step sin/cos using simple math.  Properly scaled 16 bit values
  40. >aint bad.  Test in spread sheet if you wish.  Dont forget to reset to (x=0, y=1)
  41. >after one trip around or you get a growing circle.
  42.  
  43. This is silly. Someone already posted the Bresenham algorithm which doesn't
  44. require any sines, cosines or multiplication, and approximates the circle
  45. better and better with greater resolution.
  46.  
  47. The underpinnings behind are this. You represent the circle as an implicit
  48. equation:
  49.  
  50.     f(x,y)  =  x*x  +  y*y  -  r*r  =  0
  51.  
  52. You note that if you insert an arbitrary (x,y) pair into the functin f(x,y), if
  53. you insert the coordinates of a point that is outside of the circle, you get a
  54. negative result. Inserting a point that is inside the circle, you get a
  55. postive f(x,y). The circle is the locus of points for which f(x,y) = 0.
  56.  
  57. With this in mind, you can start at the point (0, r) which you know is the
  58. point where the circle crosses the y axis and move left. You have a choice of
  59. two pixels. The one immediately to the right, (x+1,y), or the one to the right
  60. and below (x+1,y+1). You also know that your current value of f(x,y)---let us
  61. call it D---is zero (you are on the circle).
  62.  
  63. To decide which pixel to step to, you look at D. If it is negative, you know
  64. you are out too far, so you choose the pixel (x+1,y+1), else you choose the one
  65. that is (x+1,y). You plot the pixel and update your (x,y) coordinates with
  66. simple increments as decided by the value of D.
  67.  
  68. Now it turns out that you can update the value of D _precisely_ using a
  69. difference formula, without evaluating x*x + y*y - r*r, even though doing so
  70. would still be a lot better than using sines and cosine tables and such.
  71.  
  72. In this case, you increasd only the x coordinate. So the amount to be added to
  73. D to correct it is:
  74.  
  75.     D_new - D_old    =     ((x+1)^2 + y*y - r*r) - (x*x + y*y - r*r)
  76.             =    x^2 + 2x + 1 - x^2
  77.             =    2x + 1
  78.  
  79. In other words, in the C language you do a "  D += 2x + 1; ", which the
  80. compiler will quite likey strength-reduce into a shift and an add.
  81. There is no loss of accuracy unless x is so large that a shift left overflows.
  82. D never strays _too_ far from zero.
  83.  
  84. Incidentally, the steps for updating in the case that you go to the other
  85. pixel are slightly different. You must compute the change in D for both the x
  86. and y increment. and you get something like
  87.  
  88.     D_new - D_old    =    2x + 1 - 2y + 1
  89.             =    2x - 2y + 2
  90.             =     2(x - y + 1)
  91.  
  92. So, the algorithm is:
  93.  
  94. procedure circle(variable r)
  95. let x = 0
  96. let y = r
  97. let d = 0
  98.  
  99. while (x <= y)        # work over 45 degree slice
  100.     plot (x,y) pixel
  101.     plot all symmetric pixels (y,x), (-y,x), etc to get full circle
  102.     if (D < 0)    # we are out too far
  103.         x = x + 1
  104.         y = y + 1
  105.         D = D + 2 * (x - y + 1)        # optimizes, can be re-written
  106.     else
  107.         x = x + 1
  108.         D = D + 2 * x - 1
  109.     endif
  110. endwhile
  111.  
  112.  
  113. No sines, no cosines, no lookup tables, no fixed-point jumble. Lightning fast.
  114.  
  115. Of course, this algorithm quantizes the circle. It is intended for raster
  116. devices. A floating point method involving stepping through a parametric
  117. equation might me appropriate for drawing a circle on plotter-like devices or
  118. vector displays. 
  119.  
  120. There are adaptations of this algorithm with anti-aliasing, etc.
  121.  
  122. The method can be adapted for curves other than circles. I have done hyperbolas
  123. with it.
  124.  
  125. The related method of forward differences for evalating cubic functions was
  126. used with Charles Babbage's difference engine.  Today they use it to
  127. interpolate digital signals before reconstructing and sell it to you as
  128. "oversampling".
  129. -- 
  130.  
  131.